clear; dbstop if error; close all;
result_folder = [fullfile(pwd, '..') filesep];

pfid = fopen([result_folder 'numbers_reported_in_text.txt'], 'w+');

%% Tables 2 and 3 (summary statistics)
load([result_folder 'est_demand_variables_for_est.mat'], ...
    'geomkt','yr','craftdummy','quantity','price','dummy_light','dummy_lager','dummy_ale');
load([result_folder 'demand_mc_results'], ...
    'owner_name2','brand_id2');

% quantity, price, number of firms, number of products
yr_list = 2010 : 2016;

q_allmkt = sum(quantity(ismember(yr, yr_list)));

geomkt_common = unique(geomkt(yr==yr_list(1)));
for yy = 1 : length(yr_list)
    geomkt_common = intersect(geomkt_common, unique(geomkt(yr==yr_list(yy))));
end

ind_common = ismember(geomkt, geomkt_common);

% quantity
q = nan(length(yr_list), 1);
cq = nan(length(yr_list), 1);

% price
p = nan(length(yr_list), 1);
cp = nan(length(yr_list), 1);

% number of firms
N = nan(length(yr_list), 1);
cN = nan(length(yr_list), 1);

% number of products
J = nan(length(yr_list), 1);
cJ = nan(length(yr_list), 1);

% flavor share (light, lager, and ale)
F = nan(length(yr_list), 3);
cF = nan(length(yr_list), 3);

F_count = nan(length(yr_list), 3);
cF_count = nan(length(yr_list), 3);

for ky = 1 : length(yr_list)
    ind_y = yr==yr_list(ky) & ind_common;    
    ind_cy = yr==yr_list(ky) & craftdummy & ind_common;
    
    q(ky) = sum(quantity(ind_y));    
    cq(ky) = sum(quantity(ind_cy));
    
    p(ky) = sum(price(ind_y).*quantity(ind_y))/sum(quantity(ind_y));    
    cp(ky) = sum(price(ind_cy).*quantity(ind_cy))/sum(quantity(ind_cy));
    
    N(ky) = length(unique(owner_name2(ind_y)));    
    cN(ky) = length(unique(owner_name2(ind_cy)));
    
    J(ky) = length(unique(brand_id2(ind_y)));    
    cJ(ky) = length(unique(brand_id2(ind_cy)));
    
    F(ky, 1) = sum(quantity(ind_y & dummy_light))/q(ky);    
    F(ky, 2) = sum(quantity(ind_y & dummy_lager))/q(ky);    
    F(ky, 3) = sum(quantity(ind_y & dummy_ale))/q(ky);
    
    cF(ky, 1) = sum(quantity(ind_cy & dummy_light))/cq(ky);    
    cF(ky, 2) = sum(quantity(ind_cy & dummy_lager))/cq(ky);    
    cF(ky, 3) = sum(quantity(ind_cy & dummy_ale))/cq(ky);
    
    F_count(ky, 1) = length(unique(brand_id2(ind_y & dummy_light)))/J(ky);    
    F_count(ky, 2) = length(unique(brand_id2(ind_y & dummy_lager)))/J(ky);    
    F_count(ky, 3) = length(unique(brand_id2(ind_y & dummy_ale)))/J(ky);
    
    cF_count(ky, 1) = length(unique(brand_id2(ind_cy & dummy_light)))/cJ(ky);    
    cF_count(ky, 2) = length(unique(brand_id2(ind_cy & dummy_lager)))/cJ(ky);    
    cF_count(ky, 3) = length(unique(brand_id2(ind_cy & dummy_ale)))/cJ(ky);    
end


fprintf(pfid, ['Table 2 reports summary statistics based on %1.0f markets ' ...
    'present in the data every year from 2010 to 2016. '...
    'These markets account for %1.0f%% of the total quantity from all markets and years.\n\n'], ...
    length(unique(geomkt_common)), sum(q)/q_allmkt*100);

% summary of q, p, N, J (Table 2)
fileID_sum = fopen([result_folder, 'sum_stats.txt'],'w+');

fprintf(fileID_sum, '%s\t\t%s\n', 'quantity', 'annual average');
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'craft', mean(cq));
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'all', mean(q));

fprintf(fileID_sum, '%s\n', 'price');
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'craft', mean(cp));
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'all', mean(p));

fprintf(fileID_sum, '%s\n', 'number of firms');
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'craft', mean(cN));
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'all', mean(N));

fprintf(fileID_sum, '%s\n', 'number of products');
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'craft', mean(cJ));
fprintf(fileID_sum, '\t%s\t%1.0f\n', 'all', mean(J));

fclose(fileID_sum);

% summary of flavor q and J (Table 3)
fileID_flavor = fopen([result_folder, 'sum_stats_flavor.txt'],'w+');

fprintf(fileID_flavor, '%s\n', 'shares of quantities');
fprintf(fileID_flavor, '%s\t%s\t%s\t%s\n', 'share', 'light', 'lager', 'ale');
fprintf(fileID_flavor, '%s\t%1.2f%%\t%1.2f%%\t%1.2f%%\n', 'craft', mean(cF(:, 1))*100, mean(cF(:, 2))*100, mean(cF(:, 3))*100);
fprintf(fileID_flavor, '%s\t%1.2f%%\t%1.2f%%\t%1.2f%%\n', 'all', mean(F(:, 1))*100, mean(F(:, 2))*100, mean(F(:, 3))*100);

fprintf(fileID_flavor, '%s\t\t\t\t\n', 'shares of the number of products');
fprintf(fileID_flavor, '%s\t%1.2f%%\t%1.2f%%\t%1.2f%%\n', 'craft', mean(cF_count(:, 1))*100, mean(cF_count(:, 2))*100, mean(cF_count(:, 3))*100);
fprintf(fileID_flavor, '%s\t%1.2f%%\t%1.2f%%\t%1.2f%%\n', 'all', mean(F_count(:, 1))*100, mean(F_count(:, 2))*100, mean(F_count(:, 3))*100);

fclose(fileID_flavor);

%% Figure 1 (Price of Barley)
clearvars -except result_folder test_code pfid
load([result_folder 'est_demand_variables_for_est.mat'], ...
     'barley');

inflation_tbl = [2010, 1.1;
                2011, 1.07;
                2012, 1.05;
                2013, 1.03;
                2014, 1.01;
                2015, 1.01;
                2016, 1];

ym_list_all = ((2010:2016)'*100 + (1:12))';
ym_list_all = ym_list_all(:);
barley_price = nan(size(ym_list_all, 1), 1);
for k = 1 : length(ym_list_all)    
    mk = mod(ym_list_all(k), 100);    
    yk = floor((ym_list_all(k)-mk)/100);    
    barley_price(k) = barley.Value(year(barley.Date)==yk & month(barley.Date)==mk)...
        *inflation_tbl(inflation_tbl(:,1) == yk, 2);    
end

figure;
plot(barley_price);
ylabel('Global Barley Prices (2016 $/ton)', 'fontsize', 15);
xlabel('Month', 'fontsize', 15);
axis([1 84 100 350])
xticks([1 13 25 37 49 61 73 84])
xticklabels({'01/2010', '01/2011', '01/2012', '01/2013', '01/2014', '01/2015', '01/2016', '12/2016'})

print([result_folder 'barley_avg_price'], '-dpdf');

%% Tables 4, 5, 6 and Fig 3 (demand estimates, match of micro moments, elasticities, diversion ratios, and markups)
clearvars -except result_folder test_code pfid

load([result_folder 'demand_mc_results'], ...
    'para_nonlinear_transformed', 'para_linear', 'se_key_demand',...
    'logincnorm', 'elas_table', 'brand_id_by_mkt', 'brand_id2', ...
    'diver_craft', 'diver_outside',  ...
    'mkt_weight', 'mc', ...
    'diver_lager', 'diver_ale', 'diff_micro_m', 'micro_m_data');
load([result_folder 'est_demand_variables_for_est.mat'], ...
    'quantity','price', 'craftdummy', 'mktsize', 'dummy_ale','dummy_lager', 'yr', 'brand_descr');

fprintf(pfid, 'The 25%% and 75%% quantiles of the market sizes are %1.0f and %1.0f\n', quantile(mkt_weight, [0.25 0.75]));

% table 4 
para_name = {'lager', 'light', 'craft', 'import', 'ale', 'intercept', ...
        'light-craft', 'income-price', 'income-intercept', 'income-craft', 'price'};
para_nonlinear_transformed(8:10) = para_nonlinear_transformed(8:10)/logincnorm;
para_output = [para_nonlinear_transformed, para_linear(1)];
stderr_output = se_key_demand;

para_name_print = {'intercept', 'lager', 'light', 'ale', 'craft', 'light-craft', ...
        'import',  'income-intercept', 'income-craft', 'income-price', 'price'};
    
fileID_demand = fopen([result_folder, 'demand_est.txt'],'w+');
for np = 1 : length(para_name_print)       
    fprintf(fileID_demand, '%s\t', para_name_print{np});    
    ind_np = ismember(para_name, para_name_print{np});    
    fprintf(fileID_demand, '%1.2f\n\t', para_output(ind_np));    
    fprintf(fileID_demand, '(%1.2f)\n', stderr_output(ind_np));    
end
fclose(fileID_demand);

% table 5
moment_name = {'q', 'lager', 'light', 'import', 'ale', 'lager by craft',...
    'light by craft', 'craft', 'import by craft', 'ale by craft', ...
    'quantity 50-100 diff', 'quantity 100+ diff', ...
    'Pr(income<50k | craft) - mean(Pr(income<50k | other))', ...
    'Pr(income<100k | craft) - mean(Pr(income<100k | other))', ...
    'price 50-100 diff','price 100+ diff'};

moments_compare = [micro_m_data micro_m_data+diff_micro_m];
moment_name_print = {'light', 'lager', 'ale', 'import', 'craft', 'light by craft', 'lager by craft', 'q'};
n_moment_print = length(moment_name_print);

moments_compare_print = nan(n_moment_print, 2);
for k = 1 : n_moment_print   
    moments_compare_print(k, :) = moments_compare(ismember(moment_name, moment_name_print{k}), :);
end

fileID_moment = fopen([result_folder, 'demand_micromoment_est.txt'],'w+');    
for np = 1 : length(moment_name_print)       
    fprintf(fileID_moment, '%s\t', moment_name_print{np});            
    fprintf(fileID_moment, '%1.2f\t%1.2f\n', moments_compare_print(np, :));    
end
fclose(fileID_moment);

% table 6. print the elasiticities
N_top = 3; % number of the top craft products we consider
[all_brands, ~, all_ind] = unique(brand_descr);
ttq = grpstats(quantity, all_ind, 'sum');
[~, ind_q] = sort(-ttq);
top_brand = all_brands(ind_q(1:N_top));
[~, ind_choose_top] = intersect(brand_descr, top_brand);
top_id = brand_id2(ind_choose_top);

N_top_craft = 3; % number of the top craft products we consider
[craft_brands, ~, craft_ind] = unique(brand_descr(craftdummy));
ttq_craft = grpstats(quantity(craftdummy), craft_ind, 'sum');
[~, ind_craft_q] = sort(-ttq_craft);
top_craft = craft_brands(ind_craft_q(1:N_top_craft));
[~, ind_choose_top_craft] = intersect(brand_descr, top_craft);
top_craft_id = brand_id2(ind_choose_top_craft);

top_id2 = [top_craft_id; top_id];
top_elas_sum = zeros(length(top_id2), length(top_id2));
mktw_sum = zeros(length(top_id2), length(top_id2));
for k = 1 : length(elas_table)    
    ind_k = nan(length(top_id2), 1);    
    for ne = 1 : length(top_id2)       
        ind_ne = brand_id_by_mkt{k}==top_id2(ne);
        
        if any(ind_ne)
           ind_k(ne) = find(ind_ne);             
        end        
    end    
    ind_prod = ~isnan(ind_k);    
    if any(ind_prod)        
        top_elas_sum(ind_prod, ind_prod) = top_elas_sum(ind_prod, ind_prod) + ...
            elas_table{k}(ind_k(ind_prod), ind_k(ind_prod))*mkt_weight(k);        
        mktw_sum(ind_prod, ind_prod) = mktw_sum(ind_prod, ind_prod) + ...
            mkt_weight(k);        
    end    
end

elas_top = top_elas_sum./mktw_sum;
dim_elas = size(elas_top);

fileID_elas = fopen([result_folder, 'elas', '.txt'],'w+');
for n1 = 1 : dim_elas(1)
    for n2 = 1 : dim_elas(2)
        fprintf(fileID_elas, '%1.2f\t', elas_top(n1, n2));
    end
    fprintf(fileID_elas,'\n');
end
fclose(fileID_elas);

% figure 3
[craft_id2, ~, ind_craft2] = unique(brand_id2(craftdummy));
diver_craft_id2 = nan(length(craft_id2), 1);
for k = 1 : length(craft_id2)    
    ind_k = brand_id2==craft_id2(k);    
    diver_craft_id2(k) = sum(diver_craft(ind_k).*mktsize(ind_k))/sum(mktsize(ind_k));    
end
figure; histogram(diver_craft_id2, 'Normalization', 'count');
axis([0 0.4 0 50])
xlabel('Diversion Ratio to Craft Products', 'fontsize', 15)
ylabel('Number of Craft Products', 'fontsize', 15)
print([result_folder 'diversion_to_craft'], '-dpdf');

diver_non_craft = 1 - diver_outside - diver_craft;
diver_non_craft_id2 = nan(length(craft_id2), 1);
for k = 1 : length(craft_id2)
    ind_k = brand_id2==craft_id2(k);
    diver_non_craft_id2(k) = sum(diver_non_craft(ind_k).*mktsize(ind_k))/sum(mktsize(ind_k));
end
figure; histogram(diver_non_craft_id2, 'Normalization', 'count');
xlabel('Diversion Ratio to Non-Craft Products', 'fontsize', 15)
ylabel('Number of Craft Products', 'fontsize', 15)
axis([0 0.4 0 100])
print([result_folder 'diversion_to_non_craft'], '-dpdf');

mkup = price - mc;

fprintf(pfid, 'median markup of craft beers = %1.2f\n', median(mkup(craftdummy)));
fprintf(pfid, 'mean markup of craft beers = %1.2f\n', mean(mkup(craftdummy)));

mkup_craft_id2 = nan(length(craft_id2));
for k = 1 : length(craft_id2)    
    ind_k = brand_id2==craft_id2(k);    
    mkup_craft_id2(k) = sum(mkup(ind_k).*mktsize(ind_k))/sum(mktsize(ind_k));    
end
figure; histogram(mkup_craft_id2, 'Normalization', 'count');
xlabel('Markup (2016 $)', 'fontsize', 15)
ylabel('Number of Craft Products', 'fontsize', 15)
print([result_folder 'markup'], '-dpdf');

%% Figures 4 and 5 (Data Patterns Aiding Identification)
clearvars -except result_folder test_code pfid

load([result_folder 'fc_covariates2016.mat'], ...
        'dev_profit_all', 'dev_profit_one', 'entry_sim', 'craft_sim', 'in_state_sim', 'dist_sim', 'pop_profit');


% distribution of lower/upper bounds
figure; histogram(dev_profit_all./dev_profit_one, 'Normalization', 'count');
axis([0 1 0 1200]);
xlabel('Min Change in Variable Profit/Max Change in Variable Profit', 'fontsize', 10.5);
ylabel('Number of Potential Product and Market Combinations', 'fontsize', 11)
print([result_folder, 'lower_upper_bound_ratio'], '-dpdf');
fprintf(pfid, 'median ratio = %1.2f\n', median(dev_profit_all(:)./dev_profit_one(:)));

% upper bounds and entry probabilities
profit_cutoff = [0:20:200, Inf];
entry_prob_upper = zeros(length(profit_cutoff)-1, 1);
entry_prob_lower = zeros(length(profit_cutoff)-1, 1);
for k = 1 : length(profit_cutoff)-1    
    entry_prob_upper(k) = mean(entry_sim(dev_profit_one<profit_cutoff(k+1) & dev_profit_one>=profit_cutoff(k)));
    entry_prob_lower(k) = mean(entry_sim(dev_profit_all<profit_cutoff(k+1) & dev_profit_all>=profit_cutoff(k)));
end

figure; plot(0:10, entry_prob_upper, ':', 'LineWidth', 3)
hold on; plot(0:10, entry_prob_lower, 'color', [0.8500, 0.3250, 0.0980]);
axis([0, 10 0 0.8])
ylabel('Entry Probability', 'fontsize', 15)
xlabel('Group by Change in Variable Profit', 'fontsize', 15)

legend({'Max Change in Variable Profit', 'Min Change in Variable Profit'}, 'Location', 'northwest');
print([result_folder, 'lower_upper_entry_prob'], '-dpdf');

% report the quantiles in the text
fprintf(pfid, '\t%1.2f', quantile(dev_profit_one(:), [0.1 0.25 0.5 0.75 0.9])*12); fprintf(pfid, '\n');
fprintf(pfid, '\t%1.2f', quantile(dev_profit_all(:), [0.1 0.25 0.5 0.75 0.9])*12); fprintf(pfid, '\n');

% distance, entry
ub_in_state = dev_profit_one(in_state_sim>0);
lb_in_state = dev_profit_all(in_state_sim>0);
dist_in_state = dist_sim(in_state_sim>0);
y_in_state = entry_sim(in_state_sim>0);

figure; 
histogram(dist_in_state, 20, 'normalization', 'probability'); hold on; 
histogram(dist_in_state(y_in_state>0), 20, 'normalization', 'probability', 'FaceColor', 'none'); 
legend('In-State, All', 'In-State, Entered')
ylabel('Shares', 'fontsize', 15)
xlabel('Distance (km)', 'fontsize', 15)
print([result_folder, 'dist_hist_compare_state_craft'], '-dpdf');

% log(market size) and number of products
log_mktsize = log(pop_profit);
n_prod_entered = sum(entry_sim.*craft_sim, 2);

figure;
plot(log_mktsize, n_prod_entered, '.');
ylabel('Number of Craft Products', 'fontsize', 15);
xlabel('Log(Market Size)', 'fontsize', 15);
print([result_folder, 'brand_mktsize'], '-dpdf');

%% Table 7, Table SC3, Table SC1 (fixed cost estimates)
clearvars -except result_folder test_code pfid

syf = 2016; 

para_fc_name = {'Craft', 'In State', ...
    'small FC', 'medium FC', 'large FC', ...
    'small std', 'medium std', 'large std', ...
    'small scope', 'medium scope', 'large scope', 'market std'};

fc_type = {'baseline','scope','corr'};

file_name_fc = {[result_folder 'gms_cs_baseline.mat'], ...
    [result_folder 'gms_cs_scope.mat'], ...
    [result_folder 'gms_cs_corr.mat']};

if length(fc_type) ~= length(file_name_fc); error('files missing'); end

para_fc_table_to_file = cell(length(para_fc_name), length(fc_type));

for n_spec = 1 : length(fc_type)
    fc_type_n = fc_type{n_spec};   
    
    load(file_name_fc{n_spec}, 'para_gms');
    if isequal(fc_type_n, 'baseline')
        cs = para_gms';        
        para_in_state = -cs(1, :)';        
        para_craft = -cs(2, :)';        
        para_fc0 = -cs(3:5, :)';        
        para_sig = cs(6:8, :)';     
        cs_transformed = [para_craft, para_in_state, para_fc0, para_sig];        
        tab_ind = [1 2 3 4 5 6 7 8];
    elseif isequal(fc_type_n, 'scope')
        cs = para_gms';
        para_in_state = -cs(1, :)';
        para_craft = -cs(2, :)';
        para_fc0 = -cs(3:5, :)';
        para_sig = cs(6:8, :)';       
        para_first = -cs(9:11, :)';
        cs_transformed = [para_craft, para_in_state, para_fc0, para_sig, para_first];
        tab_ind = [1 2 3 4 5 6 7 8 9 10 11];
    elseif isequal(fc_type_n, 'corr')
        cs = para_gms';          
        para_in_state = -cs(1, :)';
        para_craft = -cs(2, :)';        
        para_fc0 = -cs(3:5, :)';        
        para_sig = cs(6:8, :)'; 
        para_sig_mkt = abs(cs(9, :)');
        cs_transformed = [para_craft, para_in_state, para_fc0, para_sig, para_sig_mkt];        
        tab_ind = [1 2 3 4 5 6 7 8 12];
    end
    
    cs_projected = [min(cs_transformed, [], 1)', max(cs_transformed, [], 1)'];
    
    for k = 1 : length(para_fc_name)        
        if ismember(k, tab_ind)            
            ind_k = find(tab_ind==k);            
            para_fc_table_to_file{k, n_spec} = cs_projected(ind_k, :)*12;
        end
    end
end

fileID_fc = fopen([result_folder, 'fc_est.txt'], 'w+');

fprintf(fileID_fc, '\n');

for nfc = 1 : size(para_fc_table_to_file, 1)
    fprintf(fileID_fc, '%s\t', para_fc_name{nfc});
    
    for k = 1 : size(para_fc_table_to_file(nfc, :), 2)
        if ~isempty(para_fc_table_to_file{nfc, k})
            fprintf(fileID_fc, '[%1.2f, %1.2f]\t', para_fc_table_to_file{nfc, k}(1), para_fc_table_to_file{nfc, k}(2));
        else
            fprintf(fileID_fc, '\t');
        end
    end
    fprintf(fileID_fc, '\n');
end
fclose(fileID_fc);

%% Figures 8, SC1, SC2, Tables 8, SC2, SC4 (counterfactual results)
clearvars -except result_folder test_code pfid

load([result_folder 'est_demand_variables_for_est.mat'], ...
    'geomkt', 'craftdummy', 'yr', 'brand_descr', 'quantity','price','dummy_ale','dummy_lager');

n_mktsize_grid = 15;

sycf = 2016; 

%% (1) get result files
Files = dir(fullfile(result_folder, 'cf_results_*.mat'));

n_file = length(Files);
file_names = cell(n_file, 1);
cf_output_mktid = nan(n_file, 1);
cf_output_paraid = nan(n_file, 1);
for i = 1 : length(Files)
    file_names{i} = Files(i).name;
    tmp = regexp(file_names{i},'cf_results_(\d+)_(\d+)_','tokens');
    cf_output_mktid(i) = str2double(tmp{:}{1});
    cf_output_paraid(i) = str2double(tmp{:}{2});
end   
[tmp, sortind] = sortrows([cf_output_mktid cf_output_paraid]); %sort by market then by sampled parameters
file_names = file_names(sortind);
cf_output_mktid = tmp(:,1);
cf_output_paraid = tmp(:,2); 

n_mkt = length(unique(cf_output_mktid));
n_para = length(unique(cf_output_paraid));
if n_file ~= n_mkt*n_para; error('n_file ~= n_mkt*n_para'); end

%% (2) collect results
list_variables = {'n_owner','ind_cf','config_init_by_owner','prod_id_by_owner', 'prod_attributes_m', 'para_fc0', 'para_dist', 'para_in_state', 'para_craft', 'para_profit', 'merge_owner_id', 'owner_profit1', 'owner_profit2', 'owner_profit', 'mean_p1', 'mean_p2', 'mean_p', 'ttq1', 'ttq2', 'ttq', 'cons1', 'cons2', 'cons', 'nprod1', 'nprod', 'n_firm_sim', 'merge_owner_id_in_mkt', 'added_prod', 'dropped_prod', 'ind_cf', 'pop', 'cons_surplus_no_craft', 'cons_no_de_novo', 'added_prod_by_non_merge_entry', 'added_prod_by_non_merge'};
cf_output = cell(n_file, 1);
for nf = 1 : n_file
    load(fullfile(result_folder, file_names{nf}), list_variables{:});
    
    prod_fc = cell(n_owner, 1);
    prod_dist = cell(n_owner, 1);
    prod_in_state = cell(n_owner, 1);
    for no = 1 : n_owner
        n_prod_no = length(config_init_by_owner{no});            
        in_state_no = nan(n_prod_no, 1);            
        dist_no = nan(n_prod_no, 1);            
        craft_no = nan(n_prod_no, 1);
        
        for np = 1 : n_prod_no                
            in_state_no(np) = prod_attributes_m{prod_id_by_owner{no}(np)}.in_state;                
            dist_no(np) = prod_attributes_m{prod_id_by_owner{no}(np)}.dist*in_state_no(np);                
            craft_no(np) = prod_attributes_m{prod_id_by_owner{no}(np)}.craft;
        end
       
        fc_ns = (para_fc0 + dist_no*para_dist ...
            + in_state_no*para_in_state + craft_no*para_craft)/para_profit/1000*12;%in per-capita 1000$
                   
        prod_fc{no} = fc_ns;            
        prod_dist{no} = dist_no;            
        prod_in_state{no} = in_state_no;            
    end
    
    % compute average FC of the merged products and other products
    fc_merge = mean(cell2mat(prod_fc(merge_owner_id)));        
    fc_non_merge = mean(cell2mat(prod_fc(setdiff(1:n_owner, merge_owner_id))));        
    dist_merge = mean(cell2mat(prod_dist(merge_owner_id)));        
    instate_merge = mean(cell2mat(prod_in_state(merge_owner_id)));        
    dist_non_merge = mean(cell2mat(prod_dist(setdiff(1:n_owner, merge_owner_id))));        
    instate_non_merge = mean(cell2mat(prod_in_state(setdiff(1:n_owner, merge_owner_id))));
    instate_dist_merge = sum(dist_merge.*instate_merge)/sum(instate_merge);
    instate_dist_non_merge = sum(dist_non_merge.*instate_non_merge)/sum(instate_non_merge);
    
    dist_vec = [instate_dist_merge, instate_dist_non_merge];
    fc_vec = [fc_merge, fc_non_merge];
   
    cf_outcome.dist = dist_vec;        
    cf_outcome.fc = fc_vec;        
    cf_outcome.owner_profit = [mean(owner_profit1, 1); mean(owner_profit2, 1); mean(owner_profit, 1)];        
    
    % compute the number of new entrant firms and number of exiting firms
    ind_in_mkt = false(length(config_init_by_owner), 1);
    for k = 1 : length(config_init_by_owner)
        ind_in_mkt(k) = sum(config_init_by_owner{k})>0;
    end
    firm_id_in_mkt = find(ind_in_mkt);

    n_fc_sim = length(n_firm_sim);
    n_firm_entry = nan(n_fc_sim, 1);
    for nfc = 1 : n_fc_sim
        ind_in_mkt_cf_n = false(length(config_init_by_owner), 1);
        for k = 1 : length(config_init_by_owner)
            ind_in_mkt_cf_n(k) = sum(ind_cf{nfc}{k}) > 0;
        end
        n_firm_entry(nfc) = length(setdiff(setdiff(find(ind_in_mkt_cf_n), firm_id_in_mkt), merge_owner_id));
    end
    n_firm_exit = n_firm_entry - n_firm_sim;
    cf_outcome.n_firm_entry = mean(n_firm_entry);
    cf_outcome.n_firm_exit = mean(n_firm_exit);

    % the overall price, the craft beer prices of the merged firms, the craft beer prices
    cf_outcome.mean_p = [mean_p1; mean_p2; mean(mean_p)];
    
    cf_outcome.ttq = [ttq1; ttq2; mean(ttq)];
    
    cf_outcome.cons = [cons1, cons2, mean(cons)];
    
    cf_outcome.nprod = [nprod1; mean(nprod)];        
    
    cf_outcome.n_nonmerging_firm_change = mean(n_firm_sim);        
    cf_outcome.n_firm_merge = length(merge_owner_id_in_mkt);     

    cf_outcome.added_prod = added_prod;
    cf_outcome.dropped_prod = dropped_prod;
   
    cf_outcome.size = pop;
    
    cf_outcome.cons_no_craft = cons_surplus_no_craft*pop;        
    cf_outcome.cons_no_de_novo = cons_no_de_novo*pop;
    
    cf_outcome.added_prod_by_non_merge_entry = added_prod_by_non_merge_entry;        
    cf_outcome.added_prod_by_non_merge = added_prod_by_non_merge;
    
    cf_outcome.file_names = file_names{nf};
    
    cf_output{nf} = cf_outcome;
end

%% (3) compute merger effects 
cs_repo = nan(n_file, 1);
cs_entry = nan(n_file, 1);
cs = nan(n_file, 1);
cs_craft = nan(n_file, 1);
cs_fixed = nan(n_file, 1);
n_change = nan(n_file, 1);
craft_q = nan(n_file, 1);
mktsize = nan(n_file, 1);
craft_q_merge = nan(n_file, 1);
craft_n = nan(n_file, 1);
craft_n_merge = nan(n_file, 1);
nprod_change =  nan(n_file, 1);
nfirm_change = nan(n_file, 1);        
n_nonmerging_firm_change = nan(n_file, 1);
nfirm_entry = nan(n_file, 1);
nfirm_exit = nan(n_file, 1);
dpcraft = nan(n_file, 1);
dp2 = nan(n_file, 1);
dpall = nan(n_file, 1);
dpmerge = nan(n_file, 1);
d_own_p2 = nan(n_file, 1);
n_entry = nan(n_file, 1);
n_rival = nan(n_file, 1);    
craft_profit = nan(n_file, 1);
merge_profit = nan(n_file, 1);
merge_profit2 = nan(n_file, 1);
qall = nan(n_file, 1);
qcraft = nan(n_file, 1);
qmerge = nan(n_file, 1);      
qmerge_prop = nan(n_file, 1);    

for k = 1 : n_file            
    mktsize(k) = cf_output{k}.size;
    craft_q(k) = cf_output{k}.ttq(1,2);
    craft_q_merge(k) = cf_output{k}.ttq(1, 3);
    craft_n(k) = cf_output{k}.nprod(2, 2) - cf_output{k}.nprod(1, 2);
    craft_n_merge(k) = cf_output{k}.nprod(2, 3) - cf_output{k}.nprod(1, 3);
    cs_repo(k) = mean(cf_output{k}.cons_no_de_novo) - cf_output{k}.cons(2);
    cs_entry(k) = cf_output{k}.cons(3) - mean(cf_output{k}.cons_no_de_novo);
    cs(k) = cf_output{k}.cons(3) - cf_output{k}.cons(1);
    cs_fixed(k) = cf_output{k}.cons(2) - cf_output{k}.cons(1);
    cs_craft(k) = cf_output{k}.cons(1);
    n_change(k) = cf_output{k}.nprod(2, 1) - cf_output{k}.nprod(1, 1);
    nfirm_change(k) = -cf_output{k}.n_firm_merge+cf_output{k}.n_nonmerging_firm_change+1;  
    n_nonmerging_firm_change(k) = cf_output{k}.n_nonmerging_firm_change;
    nfirm_entry(k) = cf_output{k}.n_firm_entry;
    nfirm_exit(k) = cf_output{k}.n_firm_exit;

    nprod_change(k) = (cf_output{k}.nprod(2, 2) - cf_output{k}.nprod(1, 2));
    dpcraft(k) = cf_output{k}.mean_p(3,2)-cf_output{k}.mean_p(1,2);
    dp2(k) = cf_output{k}.mean_p(2,2)-cf_output{k}.mean_p(1,2);
    
    dpall(k) = cf_output{k}.mean_p(3,1)-cf_output{k}.mean_p(1,1);            
    dpmerge(k) = cf_output{k}.mean_p(3,3)-cf_output{k}.mean_p(1,3);
    
    qall(k) = cf_output{k}.ttq(3,1) - cf_output{k}.ttq(1,1);
    qcraft(k) = cf_output{k}.ttq(3,2) - cf_output{k}.ttq(1,2);
    qmerge(k) = cf_output{k}.ttq(3,3) - cf_output{k}.ttq(1,3);
    
    qmerge_prop(k) = cf_output{k}.ttq(1,3)./cf_output{k}.ttq(1,2);
    
    d_own_p2(k) = cf_output{k}.mean_p(2,3)-cf_output{k}.mean_p(1,3);
    n_entry(k) = mean(cf_output{k}.added_prod_by_non_merge_entry);
    n_rival(k) = mean(cf_output{k}.added_prod_by_non_merge);
    
    craft_profit(k) = cf_output{k}.owner_profit(3, 1) - cf_output{k}.owner_profit(1, 1);
    merge_profit(k) = cf_output{k}.owner_profit(3, 2) - cf_output{k}.owner_profit(1, 2);
    merge_profit2(k) = cf_output{k}.owner_profit(2, 2) - cf_output{k}.owner_profit(1, 2);
end

mktsize_measure = reshape(mktsize/1e6, [n_para, n_mkt]);        
mktsize_measure = mktsize_measure(1, :)';

craft_measure = reshape(craft_q, [n_para, n_mkt]);        
craft_measure = craft_measure(1, :)';

cs_variety = cs_entry+cs_repo;        
n_merge = nprod_change-n_entry-n_rival;

% the figures: 10 groups; max, min, average
us = reshape(mktsize, [n_para n_mkt]);
us = us(1, :);  % unique market size

[~, ind_size2] = sort(us);
[~, ind_size] = sort(ind_size2);               

n_remainder = mod(n_mkt, n_mktsize_grid);        
n_rep = (n_mkt-n_remainder)/n_mktsize_grid;

gs0 = [reshape(repmat(1:n_rep, [n_mktsize_grid, 1]), [n_mktsize_grid*n_rep 1]); (n_rep+1)*ones(n_remainder, 1)];        
gs = gs0(ind_size); % group index 

% get the average across groups
merge_var = [nprod_change n_entry n_rival n_merge cs_variety./craft_q*12*8 ...
    cs./craft_q*12*8 dpcraft ...
    avgfc(:, 1) avgfc(:, 2) dp2 cs_entry./craft_q*12*8 cs_repo./craft_q*12*8, ...
    cs_fixed./craft_q*12*8 ...
    nfirm_change n_nonmerging_firm_change nfirm_entry nfirm_exit];

n_group_var = size(merge_var, 2);
merge_var = permute(reshape(merge_var, [n_para, n_mkt, n_group_var]), [2 1 3]);

group_var = nan(n_rep+1, n_para, n_group_var);
for j = 1 : (n_rep+1)            
    group_var(j, :, :) = sum(merge_var(gs==j, :, :).*us(gs==j)')/sum(us(gs==j));
end

nprod_change_group = group_var(:,:,1);        
n_entry_group = group_var(:,:,2);        
n_rival_group = group_var(:,:,3);        
n_merge_group = group_var(:,:,4);        
cs_variety_group = group_var(:,:,5);        
cs_group = group_var(:,:,6);        
dp_group = group_var(:,:,7);        
avgfc_inc_group = group_var(:, :, 8);        
avgfc_ent_group = group_var(:, :, 9);        
dp2_group = group_var(:, :, 10);        
cs_entry_group = group_var(:, :, 11);        
cs_repo_group = group_var(:, :, 12);        
cs_fixed_group = group_var(:, :, 13);
nfirm_change_group = group_var(:,:, 14);
n_nonmerging_firm_change_group = group_var(:,:,15);     
nfirm_entry_group = group_var(:,:, 16);
nfirm_exit_group = group_var(:,:, 17);

whiskerplot(nfirm_entry_group')    
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Entrants)', 'fontsize', 15)
print([result_folder, 'nfirm_entry'], '-dpdf');

whiskerplot(nfirm_exit_group')    
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Exits)', 'fontsize', 15)
ylim([0 0.5])
print([result_folder, 'nfirm_exit'], '-dpdf');

whiskerplot(n_entry_group')
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Products by Entrants)', 'fontsize', 15)
print([result_folder, 'n_entry'], '-dpdf');

whiskerplot(n_merge_group')        
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Products by Merging Firms)', 'fontsize', 15)  
print([result_folder, 'n_merge'], '-dpdf');

whiskerplot(n_rival_group')
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Products by Rival Incumbents)', 'fontsize', 15)
print([result_folder, 'n_rival'], '-dpdf');

whiskerplot(nprod_change_group')    
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('\Delta(# Products)', 'fontsize', 15)
print([result_folder, 'nprod'], '-dpdf');

whiskerplot(dp_group')
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('Price Change ($)', 'fontsize', 15)    
print([result_folder, 'price'], '-dpdf');

figure; whiskerplot(dp2_group');
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('Price Change, Variety Fixed ($)', 'fontsize', 15)
print([result_folder, 'price_change_variety_fixed'], '-dpdf');    

whiskerplot(cs_fixed_group')
xlabel('Group Sorted by Market Size', 'fontsize', 15)
ylabel('Change in CS/Pre-Merger Craft Sales ($)', 'fontsize', 15)
print([result_folder, 'average_cs_variety_fixed'], '-dpdf');

n_entry0 = n_entry;
n_rival0 = n_rival;    
nprod_change0 = nprod_change;

% counterfactual tables 
table_var_avg = [nfirm_change n_nonmerging_firm_change nprod_change0 n_merge0 ...
    n_rival0 n_entry0 dpall dpcraft dpmerge];

table_var_sum = [qall, qcraft, qmerge, cs, craft_profit, merge_profit, craft_profit+cs, ...
    cs_variety cs_entry cs_repo cs_fixed]*12/1e3;

n_group_var_avg = size(table_var_avg, 2);
table_var_avg = permute(reshape(table_var_avg, [n_para, n_mkt n_group_var_avg]), [1 3 2]);
table_group_var_avg = sum(table_var_avg.*reshape(us, [1 1 n_mkt]), 3, 'omitnan')/sum(us);        

n_group_var_sum = size(table_var_sum, 2);
table_var_sum = permute(reshape(table_var_sum, [n_para, n_mkt n_group_var_sum]), [1 3 2]);
table_group_var_sum = sum(table_var_sum, 3);

fileID_cf_change = fopen([result_folder, num2str(sycf), '_agg_outcome.txt'], 'w+');
fprintf(fileID_cf_change, '\tchange range');
fprintf(fileID_cf_change, '\n');

% average
% 1. firms
fprintf(fileID_cf_change, 'number of firms');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 1)), max(table_group_var_avg(:, 1)));
        
% 2. new entrants
fprintf(fileID_cf_change, 'number of new entrants');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 2)), max(table_group_var_avg(:, 2)));

% 3. number of products
fprintf(fileID_cf_change, 'number of products');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 3)), max(table_group_var_avg(:, 3)));

% 4. number of merging firms products
fprintf(fileID_cf_change, 'number of merging firms products');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 4)), max(table_group_var_avg(:, 4)));

% 5. rival incumbent products
fprintf(fileID_cf_change, 'number of rival firms products');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 5)), max(table_group_var_avg(:, 5)));

% 6. entrant products
fprintf(fileID_cf_change, 'number of entrant products');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 6)), max(table_group_var_avg(:, 6)));

% 7. average prices
fprintf(fileID_cf_change, 'average prices');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 7)), max(table_group_var_avg(:, 7)));

% 8. craft product prices
fprintf(fileID_cf_change, 'average craft prices');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 8)), max(table_group_var_avg(:, 8)));

% 9. merging firms craft product prices
fprintf(fileID_cf_change, 'average merging firms prices');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_avg(:, 9)), max(table_group_var_avg(:, 9)));

% 10. total quantity
fprintf(fileID_cf_change, 'total q');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 1)), max(table_group_var_sum(:, 1)));

% 11. craft quantity
fprintf(fileID_cf_change, 'total craft q');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 2)), max(table_group_var_sum(:, 2)));

% 12. craft merging quantity
fprintf(fileID_cf_change, 'total merge q');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 3)), max(table_group_var_sum(:, 3)));

% 13. consumer surplus
fprintf(fileID_cf_change, 'CS');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 4)), max(table_group_var_sum(:, 4)));

% 14. craft profits
fprintf(fileID_cf_change, 'profits');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 5)), max(table_group_var_sum(:, 5)));
        
% 15. craft merging firm profits
fprintf(fileID_cf_change, 'merging');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 6)), max(table_group_var_sum(:, 6)));

% 16. total surplus
fprintf(fileID_cf_change, 'surplus');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 7)), max(table_group_var_sum(:, 7)));

% 17. CS due to variety
fprintf(fileID_cf_change, 'CS variety');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 8)), max(table_group_var_sum(:, 8)));

% 18. CS due to entry
fprintf(fileID_cf_change, 'CS entry');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 9)), max(table_group_var_sum(:, 9)));

% 19. CS due to incumbent product adjustment
fprintf(fileID_cf_change, 'CS repo');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 10)), max(table_group_var_sum(:, 10)));

% 20. CS due to price changes
fprintf(fileID_cf_change, 'CS fixed products');
fprintf(fileID_cf_change, '\t[%1.2f, %1.2f]\n', min(table_group_var_sum(:, 11)), max(table_group_var_sum(:, 11)));


fclose(fileID_cf_change);


fclose(pfid);